# TODO: Add comment
# 
# Author: yaping
###############################################################################
#parameter: work directory, file_name_prefix, input gch, hcg/wcg files, step pace to draw line, scale, smooth or not

files<-NULL
names<-NULL
for (e in commandArgs(TRUE)) {
	ta = strsplit(e,"=",fixed=TRUE)
	if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "wd"){
			wd<-ta[[1]][2]  ## directory to make plot file
			
		}
		if(ta[[1]][1] == "prefix"){
			prefix<-ta[[1]][2]
		}
		
		if(ta[[1]][1] == "step"){
			step<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "scale"){
			scale<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "axistep"){
			axistep<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "smooth"){
			smooth<-as.numeric(ta[[1]][2])  ## 0: no smooth; 1: smooth step;
		}
		if(ta[[1]][1] == "y_scale_min"){
			y_scale_min<-as.numeric(ta[[1]][2])  
		}
		if(ta[[1]][1] == "y_scale_max"){
			y_scale_max<-as.numeric(ta[[1]][2])  
		}
		if(ta[[1]][1] == "y_step"){
			y_step<-as.numeric(ta[[1]][2])  
		}
		if(ta[[1]][1] == "output"){
			output<-ta[[1]][2] 
		}
		if(!is.na(pmatch("fn", ta[[1]][1]))){
			files<-c(files, ta[[1]][2])
		}
		if(!is.na(pmatch("name", ta[[1]][1]))){
			names<-c(names, ta[[1]][2])
		}
		#if(ta[[1]][1] == "color"){
		#	line_color<-ta[[1]][2]
		#}
		
	}
}
setwd(wd)
mainTitle=prefix
#postscript(output, paper="special", height=4, width=4)
pdf(output, paper="special", height=4, width=4)
par(oma=c(1, 1, 1, 1))
par(mar=c(1, 1, 3, 1))

files<-files[order(names)]




valueGchSum<-NULL

for (fn in 1:length(files)) { 
	gch<-read.table(files[fn],sep="\t",header=F)
	axisSeqForPlot<-seq(0-scale, scale, by=axistep)
	if(smooth == 0){ ## no smooth
		axisSeq<-seq(0-scale, scale, by=step)
		dataSeq<-seq(5+((length(gch[1,])-5)/2)-scale, 5+((length(gch[1,])-5)/2)-scale+scale * 2, by=step)
	}
	if(smooth != 0){ ## smooth
		axisSeq<-seq(0-scale, scale, by=smooth)
		dataSeq<-seq(5+((length(gch[1,])-5)/2)-scale, 5+((length(gch[1,])-5)/2)-scale+scale * 2, by=smooth)
		
	}
	valueGch<-array()
	for(i in dataSeq){
		valueGch<-cbind(valueGch,mean(colMeans(gch[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
		
	}
	valueGch<-valueGch[,2:length(valueGch[1,])]
	valueGchSum<-rbind(valueGchSum,valueGch)
	

	
}
valuePlot<-colSums(valueGchSum,na.rm=T)
plot(axisSeq,valuePlot,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="black",lty=1,font=2,lwd=3)
par(new=T)
axis(1,at=axisSeqForPlot,lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)

axis(2,at=seq(y_scale_min,y_scale_max,by=y_step),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)
title(mainTitle, cex.main = 0.6, font.main= 4, col.main= "black",xlab="Distance to elements (bp)", ylab="Density")
legend("topright","MARs(Sum)",col="black",lty=1,cex=1.2,lwd=2)
abline(v=0)
dev.off()



